Approach to ergodicity in quantum wave functions 



Bruno Eckhardt^, Shmuel Fishman^'^, Jonathan Keating'*, 
Oded Agam^, Jorg Main® and Kirsten Miiller® 
FB Physik und Institut fiir Chemie und Biologic des Meeres, 
C.v. Ossietzky Universitdt, Postfach 25 03, D-26111 Oldenburg, GERMANY 
^ Department of Physics, Technion, Haifa, 32000, ISRAEL 
"^Minerva Center for Nonlinear Physics of Complex Systems, 
Techmon, Haifa, 32000, ISRAEL 
* Department of Mathematics, The University, Manchester, M13 9PL, UK 
^ Institut fiir Theoretische Physik, Ruhr-Universitdt Bochum, 44780 Bochum, GERMANY 
^ Centre de Recherches Nucleaires, Lab. Physique Theorique, E-67037 Strasbourg, ERANCE 

According to theorems of Shnirelman and followers, in the semiclassical limit the quantum wave- 
functions of classically ergodic systems tend to the microcanonical density on the energy shell. We 
here develop a semiclassical theory that relates the rate of approach to the decay of certain classical 
fluctuations. For uniformly hyperbolic systems we find that the variance of the quantum matrix 
elements is proportional to the variance of the integral of the associated classical operator over tra- 
jectory segments of length Th, and inversely proportional to T^, where Th = hp is the Heisenberg 
time, p being the mean density of states. Since for these systems the classical variance increases 
linearly with Th, the variance of the matrix elements decays like 1/Th- For non-hyperbolic systems, 
like Hamiltonians with a mixed phase space and the stadium billiard, our results predict a slower 
decay due to sticking in marginally unstable regions. Numerical computations supporting these 
conclusions are presented for the bakers map and the hydrogen atom in a magnetic field. 
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I. INTRODUCTION 

In the semiclassical limit, quantum wavefunctions (or 
their corresponding phase space counterparts, such as the 
associated Wigner functions) are supported by classically 
invariant structures [1,2]. In integrable systems they are 
concentrated on tori [3,4] and in chaotic systems they 
tend to spread over connected components of the chaotic 
region [5] . The consequences of this for matrix elements 
have been analyzed by Shnirelman [6], Zelditch [7], Colin 
de Verdiere [8,9], and others [10,11]. Roughly speaking, 
their theorems state that in the semiclassical limit the 
matrix elements of smooth operators tend to the micro- 
canonical average. Our aim here is to study the rate at 
which they do so. 

The simplest way of defining the semiclassical limit 
is to allow a variable Planck's constant h and to con- 
sider a sequence of values h„ approaching zero such that 
for each n some fixed energy Eg is an eigenvalue of 
the Schrodinger equation. Then one has that for any 
smooth classical observable A(p, q) with which a reason- 
able quantum operator A can be associated, almost all 
diagonal matrix elements (n|j4|n) approach the classical 
microcanonical average 

{A), I = J dpdq 6(Eo - H(p, q))A(p, q)/fi (1) 

as n ^ cxD, where SI is a normalization factor such that 

Since this definition of the semiclassical limit may be 
unfamiliar, and since it might on first sight appear arti- 
ficial in a world in which k is in fact a constant, we point 



out two alternatives (see e.g. [12]). For the first, consider 
Hamiltonians homogeneous in positions and momenta, 
e.g. billiards, or systems with a suitable scaling of pa- 
rameters, such as hydrogen in a magnetic field. Then it 
is possible to absorb Planck's constant into some power 
of the energy and so to map the semiclassical limit h ^ 
into the more familiar one of increasing energy or increas- 
ing quantum numbers. 

The second alternative applies to general systems with 
non-scaling Hamiltonians. One then exploits the fact 
that the density of states near any energy Eg will be semi- 
classically large, i.e., there are many states in an interval 
over which the classical mechanics does not change very 
much. Thus it is possible to expand the actions and other 
classical quantities to first order around the reference en- 
ergy Eq. If the potential is bounded for all energies, the 
density of states will increase with increasing Eg, so that 
one can imagine covering the energy axis with intervals 
of fixed size which contain increasing numbers of states, 
but in which the classical mechanics is essentially fixed. 

Common to all three approaches is the assumption that 
eigenstates can be labelled by integers n which number 
either the values of the quantized Planck's constant, or 
the scaled energies, or the actual energies, in such a way 
that the semiclassical limit corresponds to n ^ cxd. 

The restriction on operators in the Shnirelman-type 
theorems is rather weak; it includes position and mo- 
mentum operators, and smooth functions thereof, but it 
excludes projection operators since these do not have a 
smooth classical limit. More interesting is a restriction 
to 'almost air eigenstates. This is quantified in terms 
of densities d of subsets {E„^} of states, defined as the 
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quotient of the number of states in the set to the total 
number of states, 



lim —^{jii < n} 



(2) 



The Shnirelman-type theorems [6-10] hold for subsets of 
density one. Thus they still leave room for some indi- 
vidual wave functions to show relatively large deviations 
from the average and, perhaps, to be scarred in the neigh- 
borhood of short periodic orbits [13-15]. 

Several studies have verified that quantum matrix ele- 
ments, both diagonal and off-diagonal, do indeed fiuctu- 
ate around the classical limit [16-18]. Our concern here 
is with the variance of these fiuctuations as the semiclas- 
sical limit is approached. Obviously when Ti vanishes the 
fiuctuations must also vanish if Shnirelman's result is to 
be recovered. However, this decay may be rather slow, as 
mentioned by Colin de Verdiere [9] , and it may allow for 
large deviations due to scars in wavefunctions near peri- 
odic orbits [13,14]. In order to quantify these deviations 
we propose to look at the distribution of diagonal matrix 
elements, 

PN,Mda = Prob | e[a,a + da] and N < n < N + 



(3) 

as N and M both tend to infinity. Because of the scaling 
of the density of states in rf-dimensional systems, one can 
take M ~ with j3 < {d — l)/d, so that the average 
includes an increasing number of states while the under- 
lying classical mechanics is asymptotically fixed. The 
Shnirelman-type theorems then suggest that this distri- 
bution becomes narrower as one approaches the semiclas- 
sical limit N ^ CO. 

Here we derive and test numerically a semiclassical ex- 
pression for the variance of the matrix element distribu- 
tion, relating it to the variance and a correlation function 
characterising certain classical fiuctuations. Specifically, 
we focus our attention on the variance of matrix elements 
for states that are localized in one chaotic component. 
Using recent results from periodic orbit theory [19,12], we 
improve on the arguments of Feingold and Peres [16] and 
Prosen [18], who previously related fiuctuations of matrix 
elements to classical phase space averages. Our analysis 
is somewhat similar to that of Wilkinson [20,21], but we 
are able to relax substantially his assumptions concern- 
ing the energy smoothing. We go beyond these studies 
by estimating the fiuctuations of matrix elements around 
their classical averages and relating them directly to fiuc- 
tuations of corresponding classical quantities. We ignore 
states localized in regular regions, but will consider the 
effects of classical sticking near islands in the chaotic sea. 
In particular, differences in the variances of matrix ele- 
ments in hyperbolic and nonhyperbolic systems will be 
investigated. 

The outline of the paper is as follows. In the next sec- 
tion we discuss the semiclassical derivation of the vari- 



ance of matrix element distributions from several differ- 
ent points of view. Specifically, we describe an applica- 
tion of periodic orbit theory for hyperbolic systems, a 
reformulation of this approach in terms of a correlation 
function which is also be applicable to non-hyperbolic 
systems, and a derivation based on statistical properties 
of the classical motion. Some relations to random matrix 
theory and randomness assumptions for wave functions 
are discussed in section 3. In section 4, we present nu- 
merical data for the bakers map and for hydrogen in a 
magnetic field. We conclude in section 5 with a sum- 
mary, some remarks on the stadium billiard, and addi- 
tional general comments. 



II. SEMICLASSICAL MATRIX ELEMENT 
FLUCTUATIONS 

In this section we relate fiuctuations of diagonal matrix 
elements to properties of the periodic orbits and a cor- 
relation function of the corresponding classical motion. 
Our approach is explained in section II A, and for hyper- 
bolic systems the semiclassical calculations are carried 
^put in section II B using periodic orbit theory. Unfor- 
-xunately, for nonhyperbolic systems we cannot use the 
resulting expression as it stands because of a lack of un- 
derstanding of the role of the periodic orbits in this case. 
Instead, we derive a connection to a classical correla- 
tion function as an intermediate step within the original 
framework and then relate the fiuctuations to an integral 
over this (section II C). We argue that the result can 
also be used in non-hyperbolic systems, as well as for 
chaotic components of mixed systems. The semiclassical 
matrix element distribution then depends on the decay 
of correlations in the classical system. 



A. Variances of matrix elements 

On the quantum mechanical side, we consider the ma- 
trix element weighted density of states. 



(E) = J2HA\n)6(E-E„) 



(4) 



Without loss of generality, we assume that the operator 
A has been shifted by its average, so that the matrix 
elements fiuctuate around zero. Using a trick of Berry 
[22] the variance of the matrix elements, i.e. the average 
of their square, may be obtained from the square of the 
above density. 

To avoid problems with the product of delta functions, 
we use smooth approximations, e.g. Gaussians of width 



6,(E) 
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27re 



(5) 
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Crucially, the product of two Gaussians is again delta- 

1 



function like, viz.. 



(6) 



The weighted density of states with these smoothed delta 
functions will be denoted by pi'^\E). Berry's approach 
then involves the square of p"/^^ , multiplied by a factor 
proportional to the Gaussian smoothing parameter e, 

Ki^\E) = 2V^e [pi^\E)y (7) 
= y^^y^^{n\A\n){m\A\m)2^/TT£ ■ 

n m 

6,{E - E„)6,{E - Em) ■ (8) 

If we take e smaller than the mean spacing between levels 
then, assuming there are no degeneracies, only the terms 
with E„ = Em contribute. Hence, using (6) 



En) 



(9) 



from which we can estimate the variance by averaging 
over some energy interval AE. Since the total number of 
states in such a range is pAE, p being the mean density 
of states (unweighted), the variance cr^ of the quantum 
matrix element distribution is given by 



-a(e) 



E+AE ,p 



(10) 



In the following subsections, we will derive various semi- 
classical and classical expressions for this fluctuation 
measure. 

It should be noted that the choice of other functional 
forms for the smoothing of the delta functions changes 
some of the constants in the intermediate steps, but ul- 
timately affects neither the final result for A'e^-'(£'), nor 
the fact that this has a unique, well defined limit as e ^ 0. 
A Gaussian smoothing has the advantage of being well 
localized in both the energy and time domains, so that 
there is little overlap between off-diagonal contributions 
to expressions like (8). 



B. Hyperbolic systems 

We begin with the semiclassical approximation for 
derived in [24], 

p(^^){E) = Y,{n\Mn)KE-En) 



J ^A(p,ci)6(E-H(p,ci)) + 



-l-^A,..,e-^-/\ (11) 



2irh 



where p labels periodic orbits, of which the sum extends 
over positive and negative traversals, D denotes the num- 
ber of degrees of freedom. 



A(p(t),q(t))dt 



(12) 



is the integral of the observable along the pth orbit, and 
Tp and Wp are the orbit's period and weight. As explained 
in the introduction, we shall focus our attention on states 
in the neighborhood of Eg and hence have expanded the 
orbit actions around this reference energy using Sp(E') = 
Sp{Eo) + Tp{Eo){E' - Eo), with E = E' - Eq, in which 
case the weights are given by 



g8Sp(£;o)/S-8Vp7r/2 

|det(l - Mp)|i/2 



(13) 



where Mp is the (monodromy) matrix of the linearization 
perpendicular to the orbit and i>p is the Maslov index. 

As before, we assume that the average of the operator 
vanishes so that the first term in (11) drops out. We 
then have to evaluate the square of the convolution of 
the semiclassical expression with the Gaussian smoothing 
(5). This gives as a semiclassical approximation to A'e^-* 



p' p'' 



^iiT^,-T^„ )Eln^-c^Tl,l2n^ ^-c^Tl„l2n^ ^ ^^^^ 

where the sum extends over all pairs of orbits and their 
negative traversals. 

We now claim that the main contributions to this sum 
come from the diagonal terms for which p' = p" , so that 



(15) 



where the p's label individual orbits without negative 
traversals. If the system has time-reversal symmetry, 
then orbits come in pairs with the same phase and weight, 
giving rise to a symmetry factor g = 2. In systems with- 
out time-reversal symmetry, e.g. generic systems in a 
magnetic field, there is no such pairing and g = I. 

There are two ways to justify the neglect of the off- 
diagonal contributions. One source of cancellations re- 
sults from the variation in the signs of the Ap's, which 
must be present since the average {A)ci = 0. Assuming 
that the Ap's are random, uncorrelated, and also have a 
vanishing mean one can justify (15) for any e up to the 
limit set by the requirement that the Gaussians in (8) 
do not overlap, i.e. for e < l/p (see section IID). It will 
be shown below that the diagonal approximation results 
in a well defined, e-independent value for the variance, 
despite this arbitrariness. 

Another argument, which we will take advantage of 
in the non-hyperbolic case, appeals to Berry's analysis 
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of spectral statistics [22]. If A = 1, then Ap = Tp and 
(14) is directly related to the semiclassical expression for 
the form factor of the density of states. Classical and 
semiclassical sum rules as well as numerical observations 
[25] suggest that this form factor increases with increas- 
ing t for t < Th and is constant for t >> Th , where 
Th = 2'iTh'p is the Heisenberg time. As shown in [22], the 
range t << Th is well described by the diagonal contri- 
butions to (14). Since similar behaviour is also observed 
for matrix-element- weighted form factors [23], this im- 
plies that we can use the diagonal approximation up to 
times of the order Th, and so suggests a critical size for 
e of 



1 



1 



(16) 



This approximation is naturally expected to be better for 
the CUE than the GOE. 

In the first argument, e limits the periods of contribut- 
ing orbits to Tp < Ti/e, whereas in the second some 
fc is suggested which guarantees that the interferences 
between off-diagonal contributions are killed. In hyper- 
bolic systems the choice of e is not critical, but in non- 
hyerpbolic systems with slowly decaying correlations e 
determines a cut-off in the sum (15) which does affect 
the final expression. In such cases we will fix its value 
to be that given by (16), so that the sum is effectively 
truncated at the Heisenberg time Th- Clearly, only the 
order of magnitude of e is suggested by the above argu- 
ments, not the precise value, implying some variability in 
the semiclassical estimate for the matrix elements. One 
could improve on this if the correlations between actions 
and averages Ap along periodic orbits were understood 
[25,26]. 

To evaluate (15) we have to appeal to some classical 
results [27-30]. From an analysis of the classical evolu- 
tion operator one finds that the probability of returning 
to an infinitesimal tube around a periodic orbit p after a 
time T is given by 



Pp{T) = 8{Tp-T) 



T„ 



det (Mp - 1)1 
6(Tp-T)Tp\wp\\ 



(17) 
(18) 



When summed over all orbits with periods Tp G [T,T ■ 
AT] for sufficiently large T, this satisfies 



/T + AT 
J2Pp{T)dT = AT, 



(19) 



since the periodic orbits approximate the invariant mea- 
sure. The density of orbits increases like e'^'-^/T, hf being 
the topological entropy, so that the weights \wp\^ have 
to decrease on average like e"''*"^. The integrals Ap for 
orbits in this interval will fiuctuate around zero. Assum- 
ing that correlations along the orbits decay sufficiently 
rapidly, as in hyperbolic systems, the contributions to 



the integration of the observable along the orbit will fiuc- 
tuate randomly between positive and negative values, so 
that the distribution of Ap 's for orbits with periods near 
T will be Gaussian with a variance that increases linearly 
with T, 



PT{a)da = Prob {Ap £ [a, 



J , Tp near T} 



\/2'iTaT 



(20) 



Combining the sum rule (19) with the variance implied 
by (20) gives 



T<Tp<T + AT 



aT AT 



(21) 



so that upon replacing the sum over orbits in the diagonal 
approximation (15) by an integral, we find 



A'(^)(i?) = g 



2eV^ 
'iTT^h'-' Jo 



(22) 
(23) 



where g is again the symmetry factor. With the choice 
(16) for e, the integral is effectively truncated at the 
Heisenberg time. However, because of the distribution 
(20), the final result is actually independent of e, just 
as the quantum expression itself is. The semiclassical 
estimate of the variance of the matrix elements follows 
from this after averaging over some energy interval, as 
explained before (10). The result, which is the main one 
of this section, is that 
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Th 



(24) 



This shows first of all that the variance decays as the 
inverse of the Heisenberg time, and secondly, since a is 
determined from the classical trajectories, that classical 
and quantum fiuctuations are related. 

An alternative way of writing this relation is to con- 
sider not the distribution of the integral of the observable 
along the orbits, the Ap's, but the distribution of the av- 
erages Up = Ap/Tp. This is again a Gaussian, but now 
with a variance (t'^(T) = a/T. The implication is that 
up to a symmetry related factor the widths of the quan- 
tum matrix element distribution and the distribution of 
classical averages from trajectory segments of length Th 
are the same. And since the latter decays like 1/Th, the 
quantum matrix elements narrow around the classical av- 
erage at the same rate. 



C. Correlation functions and nonhyperbolic systems 

The situation is more complicated if a system is not 
nicely hyperbolic but rather has a mixed phase space 
or marginally stable orbits. Both around islands (due 
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to trapping in the nested cantorus structure) and near 
marginally stable orbits (because of the slow escape) one 
finds increased staying times, resulting in more slowly de- 
caying correlations and possible non-Gaussian distribu- 
tions [31,32] with more slowly decaying variances. Fur- 
thermore, in the case of marginally stable orbits, the 
semiclassical weights diverge and the Gutzwiller trace 
formula used above has to be improved [33-35]. To inves- 
tigate the behavior of the connected chaotic component 
in such a system it will be useful to obtain an expression 
for the variance of matrix elements that does not depend 
explicitly on the weights of periodic orbits in any non- 
hyperbolic regions. We now derive such an expression in 
terms of a classical auto-correlation function. 

We begin by relating Ap to a correlation function along 
periodic orbits. Abbreviating the phase space argument 
of the classical observable j4(p(t), q(t)) by z(t), with an 
index indicating the periodic orbit, we have that 



At 



dti 



dr' 



Tp/2 



dt2 A(Zp(ti))A(Zp(t2)) 



Tp/2 



Tp/2 



dTA{Zp{T' + T/2))A{Zp{T'-T/2)) 



Tp/2 



drCpir) 



(25) 



where periodicity of the integrand has been exploited and 



Cp(r) 



1 



p Jo 



dr' A{zp{T' + t/2)A{zp{t' - t/2) (26) 



E 

T<T„<T + AT 



T/2 T<T„<T + AT 



T/2 



where 



Ut) = 1 - erf(2re/;j) 



(30) 



with erf(x) = ex'p(—z^)dz. The precise form of 

refiects the fact that we used a Gaussian smoothing; had 
we worked instead with Lorentzians, we would have ob- 
tained /^(t) = exp(— 2re/7}), and hence a Laplace Trans- 
form in (29). The above expression, together with the 
choice (16) for e seems to be as far as one can generally 
go in non-hyperbolic cases. It relates the quantum fiuc- 
tuations to an integral of a classical correlation function. 

The results of the previous section can be recovered 
if the correlations decay sufficiently rapidly. Then the 
integral over the correlation function tends to a constant 



lim 

T^oo 



rfrC(r) 



rfrC(r) , 



(31) 



and so in the limit e 
the estimate 



0, in which 



1, we arrive at 



(32) 



To connect this rj (eq. 31) with the variance aT of the 
Ap's (20), note that for a system with correlations decay- 
ing exponentially on a time scale A, one can write 



C(r) 



(33) 



is the auto-correlation function along the periodic orbit. 
When substituted in the orbit sum (15) one can again 
use the fact that the weighted periodic orbits approxi- 
mate the invariant density so that the average of Cp(r) 
over all orbits becomes the classical correlation function 
(for more details and a quantitative comparison in the 
hyperbolic case, see [36]). Thus we can write, in analogy 
to (21), 

T/2 

dT > ' Cp(T)Tp\ Wp\ 



Substitution into (31) thus gives rj = a. If the correla- 
tions decay more slowly, one expects rj > a. 

The relation between the variance of the matrix ele- 
ment distribution and the classical correlation function 
is now also applicable in cases of slowly decaying cor- 
relations and nonhyperbolic systems, since the decay of 
correlations is directly related to trapping in regions in 
which the periodic orbit representation of the invariant 
measure becomes suspect. If the long time behavior is 
dominated by regions where the probability to be trapped 
! for a time r decays as t~'^ , then the correlation function 
also decays in the same way: 



T/2 



rfrC(r) AT 



(27) 



C{t) ~ r- 



(34) 



where 



C{t) = {A(z(0))A(z(r))) 



(28) 

is the average correlation function as determined from 
non-periodic ergodic trajectories or, equivalently, by av- 
eraging over the invariant measure. 

The semiclassical expression (15) for the variance may 
thus be written 

-TI2 



This results from the fact that A takes similar values 
(different from the vanishing mean) in the region where 
the motion is trapped for a long time. 

The above ansatz (34) for the asymptotic behaviour of 
the correlation function is integrable for exponents 7 > 
1, leading to the same scaling for the variance of the 
matrix elements as in the hyperbolic case. For 7 < 1, 
the integral diverges and so is dominated by the effective 
cut-off at the Heisenberg time. Consequently, one has for 
the asymptotic behaviour of the semiclassical variance 



1 /"^ 
g— drC{r)f,{r) 



(29) 




T^^ for 7 > 1 
In Th /Th for 7=1 
< 1 



-7 



H 



for 7 



(35) 
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The decay of correlations has been studied for a vari- 
ety of chaotic systems. In particular, for trapping near 
elliptic islands embedded in a chaotic sea the correlations 
have been found to decay agebraically, with the expo- 
nent 7 in the range 7 ~ 1.2 — 1.5 [37-41]. This decay is 
determined by the winding numbers and cantori of the 
surrounding islands and thus may be visible only after 
exceedingly long times [41,42]. The stadium billiard falls 
into the middle category, since there 7=1 [43]. 



D. A derivation based on classical randomness 

For hyperbolic systems, correlations between trajecto- 
ries decay sufficiently fast that for long orbits different Ap 
can be considered independent random variables, leading 
to 



{{Ap,Apu)) 



dri 



dT2 {A{Zp-{Tl))A{Zpn{T2))) 



(4) 



aT 6„ 



(36) 



where ((•••)) is an average over all pairs of periodic orbits 
of period Tp G [T,T + dT] and {Ap) is the variance of the 
periodic orbit integrals as calculated from the distribu- 
tion (20). This result should also hold if the trapping is 
not too strong, since then the times at which different or- 
bits enter the trap are uncorrelated. Averaging (14) over 
periodic orbits in a small interval of periods in the vicinity 
of T, perhaps supplemented by an average over a small 
energy interval as well, leads to (15). Importantly, under 
these assumptions the diagonal approximation is not lim- 
ited to orbits of period shorter than the Heisenberg time. 
(It is worth pointing out that the difference between (14) 
and the corresponding expression for the spectral form 
factor is that the periodic orbit contributions are in this 
case proportional to Ap , and the additional randomiza- 
tion thus introduced is enough to kill the off-diagonal 
terms for all sufficiently large T.) The calculation lead- 
ing from (15) to (24) can then be repeated without the 
restriction on the periods of the orbits, resulting in the 
semiclassical matrix element variance 
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1 



dt — e 

T 



(37) 



If the correlations decay sufficiently rapidly, the distribu- 
tion of the Ap's is Gaussian and given by (20), so that 



we again arrive at 



9: 



7] 



9- 



7] 



2'irhp " Th 
in agreement with the estimate (32). 



(38) 



III. RESULTS OF THEORIES ASSUMING 
RANDOMNESS 

In this section the predictions for the variance of the 
diagonal matrix elements in the framework of theories 



that assume true randomness will be summarized. These 
will be compared with the predictions of periodic orbit 
theory derived in Sect. 2. 



A. Random Matrix Theory 

There is much evidence that many quantum proper- 
ties of chaotic systems are consistent with random ma- 
trix theory [44] , although there is no rigorous proof for 
this. It is therefore of interest to relate our predictions to 
the corresponding random matrix results. For a typical 
observable A, random matrix theory predicts [45,46] that 



(39) 



where cr\ j is the variance of the off-diagonal matrix 
elements, while g depends on the symmetry of the model 
and takes the values 1 and 2 for the GUE and GOE 
respectively. In the semiclassical limit the variance of 
the off-diagonal matrix elements is related to the classical 
correlation function by [16], 



A.off 



_q_ 

Th 



2 



dtC{t) 



(40) 



Therefore, eqs. (39) and (40) agree with eqs. (32) and (24) 
which were obtained directly from periodic orbit theory 
under the assumption that the classical correlations de- 
cay sufficiently rapidly. 



B. Random Wave Function Theory 

This theory assumes that eigenstates of chaotic sys- 
tems are Gaussian random functions [47,48]. In partic- 
ular, if we concentrate on two dimensional billiards, the 
wave function at each point q of the coordinate space 
satisfies the normalization. 



(41) 



where Sl^ is the area of the billiard. It is assumed that 
the system is invariant under time reversal and therefore 
the wave function can be considered real. Here (• • ■) gt 
denotes a statistical average over an ensemble of similar 
billiards, or over a small region in space around q. The 
correlation function between the wave function evaluated 
at different points is [47,48], 



('/'(qi)'/'(q2))st = -^Jo (fc|qi - q2| 



(42) 



where Jo is a Bessel function and k is the wave number. 
The statistical average of a diagonal matrix element re- 
duces to the microcanonical average of the observable, 
namely 
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(43) 



A. Bakers map 



which is set to vanish in the present paper. The variance 
is 

= j yrfqirfq2(t/^'(qi)t/^'(q2)).tA(qi)A(q2) 

(44) 

where we restrict ourselves to observables that depend 
only on coordinates. For Gaussian random functions sat- 
isfying (41) and (42), 

(^'(qi)^'(q2)).t = ^ (2Jo'(fc|qi - q2|) + l) (45) 

leading to the corresponding result for the variance of the 
matrix element distribution, 



_ 2_ 



j y"rfqirfq2=/o(A;|qi -q2|)A(qi)A(q2) 



In order to obtain a semiclassical expression, we use 
the (short wavelength) asymptotic form of the Bessel 

function Jo{z) ~ •y^^cos(z — 7r/4) and approximate 

cos^(z — 7r/4) by its averaged value i. The result is that 



(iqirfAq 



^(qi)^(qi + Aq) 



niT:k J J - |Aq| 

where Aq = q2 — qi • This can also be written as 



(46) 



' A,r 



2 1 
irLfj Sir 



A qQ^ qi + Aq 
rfqirfAq — — , (47) 



where Lh = Thv is the Heisenberg length, v being the 
speed. It clearly scales with the Heisenberg time in ac- 
cordance with (24) and (32). Since one can obtain (44) 
under the same semiclassical assumptions as enter our 
derivation, this result agrees with the previous expres- 
sion (O. Agam, unpublished). 



IV. NUMERICAL EXAMPLES 

Before turning to our numerical tests, we explain a 
useful classical approximation to the sum over peridoic 
orbits in (15). Since in hyperbolic systems the invari- 
ant density in phase space can be approximated by 6- 
functions on the periodic points with weights given by 
jwppTp, the above expression can, by ergodicity, also be 
computed by following a nonperiodic trajectory, which 
is then subdivided into segments of length T. This has 
been exploited in the numerical calculations below. 



As an example of a nicely hyperbolic system we take 
the bakers map, with the quantization proposed by Sara- 
ceno [49-51]. In the form given, the above formulae do 
not strictly apply to maps. Nevertheless, the time evo- 
lution is represented by a unitary operator, the traces of 
powers of which have a semiclassical periodic orbit ex- 
pansion [52-55]. In the case of the bakers map, where 
the classical phase space is bounded, the unitary matrix 
has dimension N, which corresponds to the inverse of 
Plancks constant, and all its properties can be obtained 
from traces of the first N powers. Hence, in the semiclas- 
sical expression (24), we should allow for all orbits up to 
period N. 

Because of a rapid decay of correlations in the bakers 
map, some form of the law of large numbers applies to 
averages of observables over several time steps. The vari- 
ance of the observable summed over n time steps is given 
by {Ap)(n) = an. Moreover, the distribution function 
is, to a good approximation, Gaussian, as in (20). In 
particular, for the observable A(p, x) = cos(27ra;) one has 
a = 1/2. 

Quantum eigenstates and eigenvectors were obtained 
by direct diagonalization of the unitary propagator. To 
improve on the statistics of the matrix elements, small 
groups of matrices have been collected together. If N = 
100, 200 and 400 denotes the central value, the included 
matrices are of size N , N±2 and #±4. When rescaled by 
the quantum matrix elements for cos(27ra;) also fol- 
low a Gaussian distribution, as shown in Fig. 2. Quanti- 
tatively, the variance of their distribution decreases with 
N according to 



a\ (N) 



(48) 



where a is about 1.5 . . . 1.8, the larger value correspond- 
ing to the largest N . 

The relation between the classical and quantum scaling 
with N is given by a = 4a: one factor of two coming from 
the fact that the quantum map has an antiunitary sym- 
metry and so (jf = 2, and the second factor of two because 
it also has a unitary symmetry and the two correspond- 
ing symmetry-classes contribute independently. The dif- 
ference of about 10% can probably be traced back to the 
large corrections to the semiclassical approximation in 
the bakers map due to diffraction from the discontinuity 
[55,56]. 



B. Hydrogen in a magnetic field 

As a second example we consider a smooth dynamical 
system: the hydrogen atom in a strong magnetic field 
(for reviews, see [57,58]). At the energy to be consid- 
ered here, the ergodic component covers a large fraction 
of phase space, so that no elliptic islands are visible in 
a surface of section. The tiny islands present are rarely 
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visited and only influence the very long time behaviour. 
They will not be noticable in our calculations where the 
Heisenberg time is small and fluctuations of short tra- 
jectory segments dominate. Nevertheless, the motion of 
orbits moving close to the field axis is approximately sep- 
arable and thus correlations involving these orbits may 
decay rather weakly. 

Both the classical and quantum system have scaling 
properties which ease the calculations, but which may 
obscure the relation to the previous analysis. We thus 
provide some details of the transformations involved. Af- 
ter separation of the trivial degree of freedom - a rotation 
around the magnetic field axis - the classical Hamiltonian 
for the quadratic Zeemann effect of a hydrogen atom in a 
strong magnetic field at zero angular momentum around 
the field axis becomes 



H 



2m 



1 



47r£o y/p^ + z2 8 



g ^ 2 

P 



(49) 



where (z, p) are the coordinates parallel and perpendic- 
ular to the field, (p^, pp) are the conjugate momenta, q 
is the charge of the electron, and B the magnetic field. 
Rescaling to atomic units introduces some hidden k de- 
pendencies. With the Bohr length ao = AtteoTi^ / q^m, the 
unit of momentum po = T^/o.o, the energy Eq = Ti^/maQ 
and the unit for the magnetic field Bq = niEo/qh, one 
finds the rescaled Hamiltonian (all variables denoted by 
a prime) 



H' = H/Eo 



P?+P? 



1 



+ z'2 8 



(50) 



where j = B / Bq is the dimensionless magnetic field. The 
7-dependence on the right hand side can be eliminated 
by the rescaling 



and similarly for p' , whereby 



r' = 7-2/3/' 



(51) 



Clearly one may interpret 7^/^ ^s an effective Plancks 
constant and study the semiclassical limit (as j^^^ 0) 
of matrix elements with the classical mechanics fixed. 
(It is worth pointing out that j^^^ approaches zero as 
the energy increases up to the ionization limit £" = or, 
equivalently, as the eigenvalue = y^eEg / E goes to 

infinity). This is obviously an example of the procedure 
described in the Introduction, in that one is quantizing 
(the effective) Planck's constant (7"'^''^) so that the fixed 
rescaled energy e is an eigenvalue of (54). 

The integrated density of states in this system is given 

by 



N(E,B,h) 



dpp dp;, dp dz 
2/3fi(£) 



e(E-H) 



(55) 



with 



•dp" dz" 



(27r)2 



eis-H"ip"p':,p",z"j) 



(56) 



and thus depends on energy, magnetic field and Planck's 
constant only through the dimensionless combinations 7 
and e. From this expression one can calculate the mean 
density of states as a function of energy at fixed magnetic 
field and Planck's constant, 



' dE ^0 de 



(57) 



This gives the Heisenberg time Th in original coordinates 



Th = 2TThp(E) = 2ttj 



and thus, by (53), 



2'ITj' 



-1/3^ 

de 



(58) 



(59) 



H" = 7-2/3^/^0 



P?+pf 



Vp"2 + z"2 8 



(52) 



This shows that the classical dynamics depends only on 
the rescaled energy H" = s = j~^^^E/Eo. Note that 
this quantity is independent of k. The relation between 
the original time t and the dimensionless and rescaled 
time t" in this latter system is given by 



in rescaled, dimensionless coordinates. 

For reasons of numerical convenience we choose to test 
the present theory for A = l/2r, which is smooth enough 
to allow for the application of semiclassical approxima- 
tions [24]. According to (51), the average of this variable 
over the energy shell at fixed energy and magnetic field 
scales like 72/^. Taking this scaling out, one finds a sta- 
tionary average for the matrix elements 



An = j-'''{n\A\n) 



(60) 



t = —j-H". 



(53) 



Upon quantization, eq. (52) represents a generalized 
eigenvalue problem for 7^/^ [f ^ held constant, i.e. 



1 



_(y/3)2i^_. 

^ ' 2 y/p"2 + 8 



p"' = e , 



(54) 



The classical average of the Weyl symbol of A is given by 
an integral like that for the density of states and can be 
evaluted numerically to be 



Aa = 0.259, 



(61) 



The quantum analysis is based on the calculation of 
the first 9750 positive z-parity states for e = —0.1. These 
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were obtained using the Lanczos algorithm to diagonalize 
the strongly banded Hamiltonian in a harmonic oscilla- 
tor basis. The scaled matrix elements vs. the quantized 
values of 7"^/^ are shown in Fig. 3a. The histogram, 
Fig. 3b, over all the states shows an almost Gaussian 
distribution around classical microcanonical average. 

Fig. 4 compares classical and quantum variances as 
a function of Th ■ The classical calculation is based on 
trajectory segments of length Th ■ The quantum data 
are obtained by a Gaussian smoothing of width e = 5 
in Ki^\j~^^^) (cf. eqn. 9). The classical variances de- 
cay slower than 1/Th (with an exponent of about —1/2), 
presumably due to a slower decorrelation of motion par- 
allel to the field. Thus eqns (34) and (35) have to be 
applied and the exponents found are in accord with the 
prediction. As far as the constants are concerned, the 
agreement is unexpectedly good since in this case the 
prefactors actually depend on the choice of the (arbi- 
trary) smoothing parameter e. 

V. SUMMARY AND DISCUSSION 

The main result presented here is a connection between 
the quantum fiuctuations of matrix elements around the 
classical microcanonical average and fiuctuations of clas- 
sical averages over orbit segments of length of the order 
of the Heisenberg time, as given by (24) in the hyper- 
bolic case and (29) and (10) in the nonhyperbolic case. 
These results also establish the relation (39) between di- 
agonal and off-diagonal matrix elments within periodic 
orbit theory, a direct link elusive to the approach of Peres, 
Feingold and Wilkinson [16,20]. In the case of the bak- 
ers map there is acceptable agreement between quantum 
behaviour and the semiclassical predictions, and in the 
case of the hydrogen atom in a strong magnetic field it 
is fair. The extensions proposed to non-hyperbolic sys- 
tems could not be tested in depth here for a lack of a 
sufficient number of eigenstates. The evidence from the 
hydrogen example (where the correlations did not decay 
very rapidly) suggests that the relation persists. Further 
tests should perhaps be based on maps since there one 
can probe the semiclassical limit much more deeply than 
with smooth systems. 

An investigation of non-hyperbolic systems should also 
include a discussion of the behaviour of the full distribu- 
tion of matrix elements and classical short time averages, 
rather than just the variance. The comparison of clas- 
sical and quantum variances for the quadratic Zeemann 
effect shows that even in non-ideal cases the two are of 
the same size. If there are islands of quasi-integrable 
motion, the distribution of matrix elements has a wing 
dominated by states localized on and near them. It is 
known that in such cases the classical distribution devel- 
ops algebraic tails [31,32]. Our statistics are not good 
enough to say much quantitatively about this relation. 
It would be worthwhile to see how far the connection be- 



tween the above comments about trapping in and near 
to the bouncing ball mode, or in systems with cantori, 
and fiuctuations in the matrix elements can be carried. 
This may also be of interest in mesoscopic systems, since 
one might speculate that a nice hyperbolic classical sys- 
tem will show Gaussian fiuctuations very much as in a 
disordered metal. A system with mixed phase space will 
show different fiuctuations and one might ask whether 
they are the same as those of a strongly disordered sys- 
tem showing localization. 

Along similar lines, an investigation of matrix elements 
between high lying states of the stadium billiard might be 
worthwhile, since there one expects on the classical side 
large effects due to the bouncing ball modes [59] and the 
slow decay of correlations [43] , while the quantum side is 
strongly infiuenced by scars [13,14]. 

Extensive work, in particular on the kicked rotator 
[60,61], has shown that many states are localized and 
change the distribution of matrix elements. This effect 
is clearly not captured by our result, although it is not 
obvious where the derivation has to be modified. Be- 
cause of their practical importance, for instance in low 
frequency AC conductivity, as well as noise induced dif- 
fusion (see [62,63] and references therein), this question 
deserves further attention. 

The results presented here can also be used to esti- 
mate quantitatively the error committed by a classical 
calculation, say for excitations of a molecule [64] or for 
conductivity [65]. In both cases the quantum expression 
for the correlation function will fiuctuate around the ap- 
propriate value of the classical correlation function, with 
a sigma variation given by the classical statistical varia- 
tion of trajectories up to the Heisenberg time. 

All of the calculations presented in this paper were car- 
ried out within the framework of the diagonal approxi- 
mation of the sum (15). This ignores the behaviour in 
the short-time regime, where individual orbits are impor- 
tant. Specifically, our calculations apply to the asymp- 
totic regime in which the Heisenberg time Th, the time 
scale relevant for the results (32), (35) and (35), is much 
larger than the periods of the shortest periodic orbits. 
There is, of course, a wide range of energies where the 
short orbits may be of interest, as demonstrated by the 
existence of scars [13,14]. They also infiuence thermody- 
namic properties at intermediate temperatures [66]. The 
exact contribution of these to the fiuctuations of matrix 
elements is left for further study. 
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FIG. 1. 

Scaled distributions of classical observables in the bak- 
ers map. The distributions are for n = 50, 100, 200, 400 
and 800 iterations of the map as computed from a long 
ergodic trajectory. The number of segments used to ob- 
tain the distribution is 16000 for the shortest and 1000 
for the largest iteration number. 



FIG. 2. 

Distribution of the matrix elements for the observable 
A = cos27ra; in the quantized bakers map. Five sets 
of matrix elements of size N , N ± 2 and # ± 4 have 
been superimposed and the central N values are listed. 
All matrix elements have been rescaled by \/^ and the 
histograms have been superimposed. The Gaussian has 
a width of about 1.8, close to the estimated value of a. 



FIG. 3. 

Scaled expectation values of A = l/2r vs eigenvalue 
Zn = for the lowest 9750 eigenstates in the posi- 

tive parity subspace of the quadratic Zeemann effect at 
e = —0.1 (a). The histogram in (b) shows the normal- 
ized distribution of matrix elements. Note the cluster- 
ing around the classical average of 0.259 and the almost 
Gaussian form. 



FIG. 4. 

Variance of quantum matrix elements and averages 
over classical trajectory segments vs. Heisenberg time 
Th/^tt. The full line for the quantum matrix elements 
has been obtained with a Gaussian smoothing of width 
e = 5 in eq. (9). The dashed curve was calculated from 
averages over classical trajectory segments of length T. 
The inset shows the same data on a log-log scale and 
reveals a transition from a T~^^^ dependence for short 
times to the anticipated hyperbolic law for larger 

times. 
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